##### Data analysis for media analysis with Indian CSOs ####

# Clean environment
rm(list = ls())

# Load packages
library(tidyverse)
library(readxl)
library(lfe)
library(plm)
library(stargazer)
library(car)

# @USER: PLEASE CHANGE TO YOUR WORKING DIRECTORY
setwd("C:/Users/HSmidt/Dropbox/Projekte/Paper_05 - CSO repression and Human Rights complaints/03_Code_and_Data/")
#setwd("C:/Users/chrstein/Dropbox/CSO repression and Human Rights complaints/03_Code_and_Data/")

# Load datasets
load("./10_replication_files/data_for_analysis_media_monthly.RData")
load("./10_replication_files/data_for_analysis_media_yearly.RData")

#### Define sample
table(media$Year)
media <- subset(media, Year>=2010 & Year<=2022)

table(media_monthly$Year)
media_monthly <- subset(media_monthly, Year>=2010 & Year<=2022)

#### Descriptive evidence for yearly repression variable #### 

# Show distribution of repression variable
categories <- c("No Repression", "Repression")
values <- as.data.frame(table(media$Repression))[,2]
tab <- cbind(categories,values)
tab <- as.data.frame(tab)
tab$values <- as.numeric(tab$values)
category_order <-  c("Repression", "No Repression")
tab$categories <- factor(tab$categories, levels = category_order)

##### Figure A.4 (appendix) #####
ragg::agg_png("./10_replication_files/Outputs/Figures/FigA4_AppE_hist_repression_India.png"
               , width = 297, height = 210, units = "mm", res = 600)
ggplot(tab, aes(x = categories, y = values)) +
  theme_bw() + 
  theme(plot.margin = margin(1,1,1,1, "cm")
        ,axis.text.x=element_text(size=16, angle = 0)
        ,axis.title=element_text(size=16)
        ,axis.title.y = element_text(size=16)
        ,axis.text.y = element_text(size=16),
        ,panel.border = element_blank()) + 
  geom_bar(stat = "identity", position="dodge", colour="black", alpha=0.2 ) + xlab("") + ylab("Number of CSO-years") 
dev.off()


# Create running variable of years until first appearance
media <- media %>%
  group_by(CSO) %>%
  mutate(count_year_firstcomm = ifelse(first_appearance == 1, 0, row_number() - which.max(first_appearance == 1))) %>% 
  ungroup()

# Problem: different counts may appear in different frequencies
table(media$count_year_firstcomm)
### Yes, so we need average values

# Get averages for each time point in relation to the first Communication
aggdat <- media %>% group_by(count_year_firstcomm) %>% 
  summarize(avg_repression = mean(Repression, na.rm = T))

# Plot average repression values
##pdf(file = "./outputs/media_reported_repression_India.pdf", width= 10)
ggplot(aggdat, aes(x= count_year_firstcomm, y = avg_repression)) +
  geom_col(position="dodge", colour="black", alpha=0.2 ) + geom_vline(aes(xintercept = 0), lty = 4) + xlab("Years from first UNSP communication") +
  ylab("Average value of media-reported repression against Indian CSOs") + theme_bw()  + 
  scale_x_continuous(breaks = seq(-10, 13, 2)) +
  theme_bw() + 
  theme(plot.margin = margin(1,1,1,1, "cm")
        ,axis.text.x=element_text(size=12, angle = 0)
        ,axis.title=element_text(size=12)
        ,panel.border = element_blank()) 
##dev.off()

#### Descriptive evidence for monthly repression variable #### 

# Show distribution of repression variable of monthly analysis
categories <- c("No Repression", "Repression")
values <- as.data.frame(table(media_monthly$Repression))[,2]
tab <- cbind(categories,values)
tab <- as.data.frame(tab)
tab$values <- as.numeric(tab$values)
category_order <-  c("Repression", "No Repression")
tab$categories <- factor(tab$categories, levels = category_order)

##pdf(file = "./outputs/hist_repression_monthly_India.pdf")
ggplot(tab, aes(x = categories, y = values)) +
  theme_bw() + 
  theme(plot.margin = margin(1,1,1,1, "cm")
        ,axis.text.x=element_text(size=12, angle = 0)
        ,axis.title=element_text(size=12)
        ,panel.border = element_blank()) + 
  geom_bar(stat = "identity",position="dodge", colour="black", alpha=0.2 ) + xlab("") + ylab("Number of CSO-years") 
##dev.off()

# Create running variable of years until first appearance
media_monthly <- media_monthly %>%
  group_by(CSO) %>%
  mutate(count_month_firstcomm = ifelse(First_Communication == 1, 0, row_number() - which.max(First_Communication == 1))) %>% 
  ungroup()

# Get averages for each time point in relation to the first Communication
aggdat_monthly <- media_monthly %>% group_by(count_month_firstcomm) %>% 
  summarize(avg_repression = mean(Repression, na.rm = T))

# Plot average repression values conditional on months before/since first communication
##pdf(file = "./outputs/media_reported_repression_India_monthly.pdf", width= 10)
ggplot(aggdat_monthly, aes(x= count_month_firstcomm, y = avg_repression)) +
  geom_col(position="dodge", colour="black", alpha=0.2 ) + geom_vline(aes(xintercept = 0), lty = 4) + xlab("Months from first UNSP communication") +
  ylab("Average value of media-reported repression against Indian CSOs") + theme_bw()  + 
  scale_x_continuous(breaks = seq(-10, 10, 2)) +
  theme_bw() + 
  theme(plot.margin = margin(1,1,1,1, "cm")
        ,axis.text.x=element_text(size=12, angle = 0)
        ,axis.title=element_text(size=12)
        ,panel.border = element_blank()) 
##dev.off()

# Test model for multi-collinearity
mtest <- lm(Repression ~ v2meharjrn + Homepage_English + Facebook + Twitter + human_rights_cso, data=media)
summary(mtest)
vif(mtest)
vif_values <- vif(mtest)

# Prepare vector of controls 
controls <- paste0(" + v2meharjrn + Homepage_English + Facebook + Twitter + human_rights_cso")


#### Run models for Yearly variable #### 
media$years_since <- media$treatment_years_since 

f0 <- as.formula( paste0("Repression  ~ treatment | 0 | 0 | CSO"))
f1 <- as.formula( paste0("Repression ~  treatment | CSO | 0 | CSO"))
f2 <- as.formula( paste0("Repression ~  treatment + years_since", controls,  "| 0 | 0 | CSO"))
f3 <- as.formula( paste0("Repression ~  treatment + years_since ", controls, " + lag(Repression,1) | lag(Repression,2:99)") )

m0 <- felm(formula = f0, data = media )
m1 <- felm(formula = f1, data = media )
m2 <- felm(formula = f2, data = media )
m3 <- pgmm(formula = f3, data = media, index=c("CSO", "Year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
summary(m0)
summary(m1)
summary(m2)
summary(m3)


# Reduced form of Table
custom_significance <- c("\\dag", "*", "**", "***")
star_reduced <- stargazer(m0, m1, m2, m3
                          , type="latex"
                          , title = "Media-reported CSO repression and UNSP communications"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)) ,
                                             c("CSO fixed effects", "No", "Yes", "No", "Yes"),
                                             c("Controls", "No", "No", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes"))
                          , keep =  c("treatment")
                          , covariate.labels = c("Post-communication")
                          , column.labels = c("(1)","(2)","(3)","(4)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "DV: Media-reported repression"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:mediaanalysis_yearly"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{12cm}{ \\textit{Notes:} The models in columns 1, 2, and 3 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3 and 4, we control for the harassment of journalists on the country-level, whether a CSO is primarily dedicated to human rights, whether it has an English homepage and appearances on Facebook and Twitter, and the number of years since the first communication targeting the CSO. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_reduced[grepl("Note",star_reduced)] <- note.latex
cat(star_reduced, sep = "\n")
#write(star_reduced, "./outputs/mediaanalysis_yearly_reduced.tex")


##### Table A.21 for Appendix E #####
star_full <- stargazer(m0, m1, m2, m3
                          , type="latex"
                          , title = "Media-reported CSO repression and UNSP communications"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)) ,
                                             c("CSO fixed effects", "No", "Yes", "No", "Yes"),
                                             c("Controls", "No", "No", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes")) 
                          , covariate.labels = c("Post-communication"
                                                 , "Count years post-treatment"
                                                 , "Harassment of journalists"
                                                 , "English homepage"
                                                 , "Facebook page"
                                                 , "Twitter page"
                                                 , "Human rights NGO"
                                                 , "Media-reported CSO repression (lag 1 yr)"
                                                 , "Constant")
                          , column.labels = c("(1)","(2)","(3)","(4)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "Dependent variable: Media-reported repression"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:mediaanalysis_yearly_app"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{12cm}{\\textit{Notes:} The models in columns 1, 2, and 3 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3 and 4, we control for the harassment of journalists on the country-level, whether a CSO is primarily dedicated to human rights, whether it has an English homepage and appearances on Facebook and Twitter, and the number of years since the first communication targeting the CSO. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/mediaanalysis_yearly_app.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA21_AppE_mediaanalysis_yearly.tex")




## Substantive interpretation

# Effects
effect_sizes_csos <- matrix(NA, nrow=4, ncol=5)
coef_m0 <- MASS::mvrnorm(1000, coefficients(m0), vcov(m0))[,2] / sd( model.frame(m0)[,c("Repression")] )
effect_sizes_csos[1,] <- quantile(coef_m0, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m1 <- MASS::mvrnorm(1000, coefficients(m1), vcov(m1))[,1] / sd( model.frame(m1)[,c("Repression")] )
effect_sizes_csos[2,] <- quantile(coef_m1, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m2 <- MASS::mvrnorm(1000, coefficients(m2), vcov(m2))[,2] / sd( model.frame(m2)[,c("Repression")] )
effect_sizes_csos[3,] <- quantile(coef_m2, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m3 <- MASS::mvrnorm(1000, coefficients(m3), vcov(m3))[,1] / sd( model.frame(m2)[,c("Repression")] )
effect_sizes_csos[4,] <- quantile(coef_m3, c(0.025, 0.05, 0.5, 0.95, 0.975))

# Compared to effect of armed conflict
names(model.frame(m3))
coefficients(m3)[1] / sd( model.frame(m2)[,c("Repression")] ) 
coefficients(m3)[6] / sd( model.frame(m2)[,c("Repression")] ) 


##### Figure 7 (main ms) #####
ragg::agg_png("./10_replication_files/Outputs/Figures/Fig7_mainms_marg_effect_plot_crosscso_media.png"
               , width = 297, height = 210, units = "mm", res = 600)

par(mar = c(5.1, 5, 2.1, 1.8), mgp = c(3.5, 1, 0), xpd = TRUE)
plot( y = c(1,2,3,4)
      , x = rev(effect_sizes_csos[,3])*100
      , pch= rev(c(16,15,17,18,20))
      , cex=1.5
      , type="p"
      , cex.lab = 1.5
      , axes= F, xlab = "Effect of post-UNSP communication period (as % st. dev. of media-rep. repression)"
      , ylab = ""
      , ylim = c(0.5, 4.5), xlim = c(-50,100))
arrows( x0 = rev(effect_sizes_csos[,1]*100), x1 = rev(effect_sizes_csos[,5]*100)
        , y0 = c(1,2,3,4),      y1 = c(1,2,3,4)
        , angle=90, length=0, lwd=1)
arrows( x0 = rev(effect_sizes_csos[,2]*100), x1 = rev(effect_sizes_csos[,4]*100)
        , y0 = c(1,2,3,4),      y1 = c(1,2,3,4)
        , angle=90, length=0, lwd=2)
lines(y=c(0.5,4.5), x = c(0,0), lty="dashed")

axis(1, at = seq(-50, 100, by=50), labels = seq(-50, 100, by=50))
text(x = rep(-60,5), y=c(1,2,3,4), labels=c( "With CSO FEs, \ncontrols, and instrumented \nlagged DV"
                                             , "With CSO FEs and \ncontrols"
                                             , "With CSO FEs"
                                             , "Bivariate regression"), cex=1.5, pos = 4)
text(x= -60, y=4.5, labels = "Models:", font=2, cex=1.5, pos=4)
dev.off()


#*############################################*#
#### Appendix for media analysis MONTHLY ####
#*############################################*#

controls_monthly <- c("Homepage_English", "Facebook", "Twitter", "human_rights_cso")
media_monthly %>% filter(!is.na(Repression)) %>% summarize(Homepage_English = mean(Homepage_English))
media_monthly %>% filter(!is.na(Repression)) %>% summarize(Facebook = mean(Homepage_English))

media_monthly$months_since <- media_monthly$treatment_months_since 

#### Run models for Monthly variable #### 
f0 <- as.formula( paste0("Repression  ~ treatment + months_since | 0 | 0 | CSO"))
f1 <- as.formula( paste0("Repression ~  treatment + months_since | CSO | 0 | CSO"))
f2 <- as.formula( paste0("Repression ~  treatment + months_since  ", controls,  "| 0 | 0 | CSO"))
f3 <- as.formula( paste0("Repression ~  treatment + months_since  ", controls, " + lag(Repression,1) | lag(Repression,2:99)") )

m0 <- felm(formula = f0, data = media_monthly)
m1 <- felm(formula = f1, data = media_monthly)
m2 <- felm(formula = f2, data = media_monthly)
m3 <- pgmm(formula = f3, data = media_monthly, index=c("CSO", "Month"), transformation = 'ld',  model = 'twostep', effect = 'individual')

summary(m0)
summary(m1)
summary(m2)
summary(m3)

##### Table A.22 for Appendix E #####
star_reduced <- stargazer(m0, m1, m2, m3
                          , type="latex"
                          , title = "Media-reported CSO repression and UNSP communications (monthly analysis)"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)) ,
                                             c("CSO fixed effects", "No", "Yes", "No", "Yes"),
                                             c("Controls", "No", "No", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes"))
                          , covariate.labels = c("Post-communication"
                                                 , "Count months post-treatment"
                                                 , "Harassment of journalists"
                                                 , "English homepage"
                                                 , "Facebook page"
                                                 , "Twitter page"
                                                 , "Human rights NGO"
                                                 , "Media-reported CSO repression (lag 1 yr)"
                                                 , "Constant")
                          , column.labels = c("(1)","(2)","(3)","(4)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "Dependent variable: Media-reported repression"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:mediaanalysis_monthly"
                          , table.placement="ht!"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{12cm}{ \\textit{Notes:} The models in columns 1, 2, and 3 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3 and 4, we control for the harassment of journalists on the country-level, whether a CSO is primarily dedicated to human rights and whether it has an English homepage and appearances on Facebook and Twitter. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_reduced[grepl("Note",star_reduced)] <- note.latex
cat(star_reduced, sep = "\n")
#write(star_reduced, "./outputs/mediaanalysis_monthly.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA22_AppE_mediaanalysis_monthly.tex")


#### Re-run the analysis with only domestic CSOs from India ####

media_domestic <- media %>% filter(!CSO %in% c("Amnesty International India", 
                                               "Asian Federation Against Involuntary Disappearances",
                                               "Asian Forum for Human Rights and Development",
                                               "Committee on Human Rights", "Greenpeace India",
                                               "Internet Watch Foundation", "Kurdish Red Crescent",
                                               "South Asian Community Centre for Education and Research",
                                               "UNICEF") )


#### Run models for Yearly variable #### 
media$years_since <- media$treatment_years_since 

f0 <- as.formula( paste0("Repression  ~ treatment | 0 | 0 | CSO"))
f1 <- as.formula( paste0("Repression ~  treatment | CSO | 0 | CSO"))
f2 <- as.formula( paste0("Repression ~  treatment + years_since", controls,  "| 0 | 0 | CSO"))
f3 <- as.formula( paste0("Repression ~  treatment + years_since ", controls, " + lag(Repression,1) | lag(Repression,2:99)") )

m0 <- felm(formula = f0, data = media_domestic)
m1 <- felm(formula = f1, data = media_domestic)
m2 <- felm(formula = f2, data = media_domestic )
m3 <- pgmm(formula = f3, data = media_domestic, index=c("CSO", "Year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
summary(m0)
summary(m1)
summary(m2)
summary(m3)


# Create table for paper
custom_significance <- c("\\dag", "*", "**", "***")
star_reduced <- stargazer(m0, m1, m2, m3
                          , type="latex"
                          , title = "Media-reported CSO repression and UNSP communications"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)) ,
                                             c("CSO fixed effects", "No", "Yes", "No", "Yes"),
                                             c("Controls", "No", "No", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes"))
                          , keep =  c("treatment")
                          , covariate.labels = c("Post-communication")
                          , column.labels = c("(1)","(2)","(3)","(4)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "DV: Media-reported repression"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:mediaanalysis_onlydomestic_yearly"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{12cm}{ \\textit{Notes:} The models in columns 1, 2, and 3 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3 and 4, we control for the harassment of journalists on the country-level, whether a CSO is primarily dedicated to human rights, whether it has an English homepage and appearances on Facebook and Twitter, and the number of years since the first communication targeting the CSO. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_reduced[grepl("Note",star_reduced)] <- note.latex
cat(star_reduced, sep = "\n")
#write(star_reduced, "./outputs/mediaanalysis_onlydomesticsos_yearly.tex")

# Create table for Appendix
star_full <- stargazer(m0, m1, m2, m3
                       , type="latex"
                       , title = "Media-reported CSO repression and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , notes.align ="l"
                       , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)) ,
                                          c("CSO fixed effects", "No", "Yes", "No", "Yes"),
                                          c("Controls", "No", "No", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes")) 
                       , covariate.labels = c("Post-communication"
                                              , "Count years post-treatment"
                                              , "Harassment of journalists"
                                              , "English homepage"
                                              , "Facebook page"
                                              , "Twitter page"
                                              , "Human rights NGO"
                                              , "Media-reported CSO repression (lag 1 yr)"
                                              , "Constant")
                       , column.labels = c("(1)","(2)","(3)","(4)")
                       , dep.var.labels.include = F
                       , dep.var.caption = "Dependent variable: Media-reported repression"
                       , column.sep.width = "1pt"
                       , align=F
                       , model.numbers=F, model.names=F
                       , label="tab:mediaanalysis_onlydomestic_yearly_app"
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{12cm}{\\textit{Notes:} The models in columns 1, 2, and 3 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3 and 4, we control for the harassment of journalists on the country-level, whether a CSO is primarily dedicated to human rights, whether it has an English homepage and appearances on Facebook and Twitter, and the number of years since the first communication targeting the CSO. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./outputs/mediaanalysis_onlydomestic_yearly_app.tex")

## Robustness test: How much do lead treatments and lag treatments predict repression
f0 <- as.formula( paste0("Repression  ~ dplyr::lead(treatment, 1) + treatment + dplyr::lag(treatment, 1) | 0 | 0 | CSO"))
f1 <- as.formula( paste0("Repression ~  dplyr::lead(treatment, 1) + treatment + dplyr::lag(treatment, 1)  | CSO | 0 | CSO"))
f2 <- as.formula( paste0("Repression ~  dplyr::lead(treatment, 1) + treatment + dplyr::lag(treatment, 1)  + dplyr::lead(later_appearances, 1) + dplyr::lag(later_appearances,1) | CSO | 0 | CSO"))
f3 <- as.formula( paste0("Repression ~  dplyr::lead(treatment, 1) + treatment + dplyr::lag(treatment, 1) + lag(Repression,1) | lag(Repression,2:99) ") )
m0 <- felm(formula = f0, data = media_monthly)
m1 <- felm(formula = f1, data = media_monthly)
m2 <- felm(formula = f2, data = media_monthly)
m3 <- pgmm(formula = f3, data = media_monthly, index=c("CSO", "Month"), transformation = 'ld',  model = 'twostep', effect = 'individual')
summary(m0)
summary(m1)
summary(m2)
summary(m3)

#### exit ###
